Expansion of a mesoscopic Fermi system from a harmonic trap 
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We study quantum dynamics of an atomic Fermi system with a finite number of particles, N, after it is 
released from a harmonic trapping potential. We consider two different initial states: The Fermi sea state and 
the paired state described by the projection of the grand-canonical BCS wave function to the subspace with a 
fixed number of particles. In the former case, we derive exact and simple analytic expressions for the dynamics 
of particle density and density-density correlation functions, taking into account the level quantization and 
possible anisotropy of the trap. In the latter case of a paired state, we obtain analytic expressions for the density 
and its correlators in the leading order with respect to the ratio of the trap frequency and the superconducting 
gap (the ratio assumed small). We discuss several dynamic features, such as time evolution of the peak due to 
pair correlations, which may be used to distinguish between the Fermi sea and the paired state. 
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Introduction — Measurements of density distributions and 
noise correlations in time-of-flight expansion have proven to 
be a very powerful direct probe of the quantum states of cold 
atomic systems. In bosonic systems, density distributions 
contain valuable information about the initial quantum state 
and allow one to observe directly Bose-Einstein condensates. 
In cold Fermi systems, the time-of-flight density profile of a 
Fermi liquid is expected to be qualitatively similar to the one 
of a paired state. Altman et al. fljj] have shown that to identify 
the existence of a BCS condensate, one should probe noise 
correlations, which in the condensed state would acquire a 
peak at the opposite momenta corresponding to Cooper pair- 
ing. Greiner et al. J2l have experimentally observed similar 
pairing correlations of fermionic atoms on the BEC side of a 
Feshbach resonance, showing the experimental feasibility of 
the proposed method. However, so far there have been no 
direct observations of pairing correlations in the BCS state, 
which was detected in Ref. [3] by different means. 

Most previous theoretical studies have concentrated on the 
dynamic properties of density and its correlators, which fol- 
lowed from the grand canonical wave-functions. However, 
all real experimentally studied systems are finite and non- 
uniform due to the trapping potential. Moreover, the grand 
canonical description is strictly speaking not appropriate in 
this context, since non-trivial density fluctuations may exist 
only if the system is finite and formally vanish in thermo- 
dynamic limit. It is therefore important to consider expan- 
sion dynamics of a finite system in a realistic trapping po- 
tential. We should mention here that there exist a number of 
exact results for a Fermi system in a harmonic potential: See 
e.g., Ref. 0], which considered thermodynamic properties of 
a mesoscopic Fermi gas and Refs. J^, 0, 0, B 0] where den- 
sity and energy distributions were studied. Also, Bruun and 
Heiselberg [ 10] investigated the properties of a paired state of 
trapped fermions. 

In this work we extend these studies to the dynamic regime 
of an expanding Fermi cloud. We explicitly calculate density 
distributions and noise correlations using the canonical for- 
mulation and taking into account discrete energy levels and 



a finite number of particles in a harmonic trap. In the case 
of a Fermi sea state, we derive exact and simple analytic ex- 
pressions for the density and its correlators. To describe a 
paired state, we use the grand canonical BCS wave-function 
projected onto the subspace with a fixed number of particles - 
projected BCS (PBCS) state. This wave-function was written 
explicitly in the original paper of Bardeen, Cooper, and Shri- 
effer [11] and then used in the context of nanoscale supercon- 
ductors 111 211 . We should mention here that there exists an ex- 
act solution of the reduced BCS Hamiltonian due to Richard- 
son |l3, 14], who found a set of equations, which determine 
the energy levels in a finite Fermi system with attractive inter- 
actions and constructed an exact wave-function in this case. 
However, to use this exact solution is difficult due to the com- 
plicated structure of the Richardson's equations. Moreover, it 
has been shown |12] that the exact Richardson's wave func- 
tion reduces to the PBCS state if the superconducting gap A 
is much larger than the level spacing (in fact the BCS Ansatz 
remains a good approximation even if the ratio of the gap and 
level spacing is of order unity Il2lll5l0 . Below, we consider 
the latter limit, which in the context of a trapped Fermi gas 
corresponds to A » to, where u> is the frequency of the trap- 
ping potential. 

Time-dependent density — Let us consider a Fermi system 
in a harmonic trap in a quantum many-body state |<D). We are 
interested in the dynamics of density and its correlators after 
the harmonic potential and interactions, if present, are turned 
off at t — 0. For simplicity, we start with the one-dimensional 
case. We will see that generalization of all results to higher 
dimensions is straightforward. The time-dependent density is 



(n(x,t)) = <O|# + (# + (x)^)£/ (0l<I>>, 



(1) 



where Uo(f) = exp(— 1'2 



tj is the evolution operator 



of the free system, E p = p /(2m), and ip(x) = i[r n (x)a n is 
the field operator, where iff„(x) and a„ are the oscillator wave- 
function and the annihilation operator in the harmonic poten- 
tial corresponding to the n-th level. Introducing the Fourier 
transform of the oscillator wave-functions i/t n {p), we can write 



2 



the density in the following form 

<**-S/i5^«('*!)*.('-!)« 

n,m 

(2) 

We note that if we neglect the g-dependence of the wave func- 
tions (which is a good approximation if t » or 1 ), we recover 
the "natural" result for the density [ 1 ] 



We therefore reproduce the classical formula of Ref. lfl6ll . Ac- 
cording to Eq. 0, the ratio crosses over from the initial value 
wjwi < 1 at t — to a final constant value 1 at t — > oo. We 
mention here that the dynamics of the single particle wave- 
functions corresponding to low-lying levels of the harmonic 
potential are different: The ratio of the spreads of the cor- 



responding wave-packets behaves as fy- = 

F 6 F & v^r -y/i+^F 

crosses over from the initial value y/cj z /a> ± < 1 at t — to a 
final constant value Ja>x/o> z > 1 at f — > oo. 

In the case of an isotropic three-dimensional trap, we can 
rewrite Eq. © in the form 



and 



(»>*^J(P)WP)(aX> 



P=mx/t 



We also note that Eq. (12) is exact and valid for any many -body 
state |<P), assuming that interactions are absent at t > 0. 
Fermi gas — Now let us consider a Fermi sea state at zero 

AM 

temperature, which means that |<D) = Yl ««l vac )- We find 



<»(r,f)> = J] (2/ +l)^[r, «,(?)], 



(8) 



where R„i(r,co) is the radial wave-function of a 
am _ „ , , , / > / '\ harmonic potential with a frequency a>, R„i = 

(n(x, r)> = V - - r | ^ ^ fl, ML ffJ -£= Cte-Wir/tfF [-(« - 0/2, i + 3/2, -r 2 /^ 2 ] with F(-) 
^2"n\yPhnojJ 2n 2n \^>J \ V^/ being the hypergeometric function. 

We also mention here that generalization of the obtained 
results to the case of a finite temperature is straightforward. 
Indeed, Eq. ([8]) describes a superposition of single-particle 
wave-functions corresponding to the occupied states. At zero 
temperature all states below the Fermi level are occupied with 
the equal probability. At finite temperatures, we have to 
weight each state with the probability determined by the Fermi 

distribution function f(E,T) = {l + exp [ £ ~^. (r) j| , where 
p(T) can be found from the relation 9nf(En) = N with 
q n — (n + l)(n + 2)/2 being the degeneracy of the n-th level. 
Thus, it is natural to write the finite temperature result as 



/raw/ \ y/ma> 

x e i( P-p"> x e i(E p'- E i>X' + i;) (3) 

The integrals over p and p' are table integrals and we obtain 

AT— 1 



(n(x,t)) = 2^[x,w(*)L 



(4) 



n=0 



where \\i n [x, u)(t)] is the oscillator wave function with a 
"rescaled" frequency u>(i) — cj/(1 + u) 2 t 2 ). The sum in Eq. 
can be calculated exactly and we find 



(n(x, t)) 



2 N (N- 1)! 



m,\-\-H 



L N\~\ ~ Ll N+ 1 I — I Hn- 1 | — 



(5) 

where we introduced £ = Vl + u 2 t 2 / y[mu>. We see that the 
behavior of the density is determined by the wave-functions 
near the "Fermi level" rip = N — 1. One can also check by 
explicit calculation that the integral over x of Eq. <(5j naturally 
reproduces the total number of particles in the system. It is 
also easy to generalize the result to the three-dimensional 
case. In the case of an anisotropic trap with frequencies u) x ,y,z 
in the corresponding directions, we find 

<n(r, r)> = J] ^ [x, a) x (t)] ^ [y, co y (t)] 0* [z, co z (t)] ,(6) 



where Ep is the Fermi level, E n ^„ i)M „ = a> (n x + n y + n, + 3/2), 
and (jD a (t) - co a /(l + co 2 t 2 ). This result allows us to deter- 
mine the dynamic ratio of the radii of the time-dependent den- 
sity distribution. E.g., let us assume that a> x = Uy — a> ± > a>,. 
The dynamic aspect ratio is determined by the ratio of the 
spreads of the wave-functions corresponding to the fermions 
occupying the highest possible levels = Ep/a> xz . 



<»(r, t- T)) (2/ + 1)4 [r, co(t)] f(E,„ T). 

n<np,l 



(9) 



Similarly, one can obtain the density-density noise correla- 
tions of a mesoscopic Fermi gas. Below, we just give a finite 
result of the calculation in the one-dimensional case at T — 

(n(xi,t)n(x2,t)) = S(x\ - X2){n{x\,t)) + {n(x\,t)){n{xz,t)) 

7 . 7 

n 2 



e ^ [H N ^(f)H N (f)-H N (f)H N ^(f)] 
n(N- l)\ 2 2 2N ( Xl - x 2 ) 2 



(10) 



(7) 



where {n(x x l , f)) is given by Eq. Q. We note that in thermo- 
dynamic limit the density correlations vanish. We also men- 
tion that higher-order density correlation functions of a Fermi 
gas can be obtained analytically and have a structure similar 
to Eq. CUll. 

Paired state — Now we consider density dynamics of an 
expanding finite Fermi system initially in the paired state. If 
the system were infinite, we could use the grand canonical 
BCS wave-function. However, the latter is not appropriate 
in the given context since the grand canonical BCS Ansatz 
does not conserve the number of particles. A superconducting 
wave-function with a fixed number of particles was discussed 
by Bardeen, Cooper, and Schrieffer in their seminal paper Hill 
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and later in great detail by von Delft and others in the 
context of ultrasmall superconducting grains. This issue was 



also investigated by Richardson | T3I 



who found an exact 



many body wave function of the reduced Hamiltonian 



6 = ^(5, 



(11) 



where \n) and \n T ) = f\n) are pairs of single particle states 
related to each other by time-reversal, the sum runs over all 
interacting single-particle levels, and s n are the energies of 
the levels. Model Hamiltonian ( fTTT i has been used to describe 
the properties of nuclei in the nuclear pairing model as well 
as small superconducting grains with a finite number of parti- 
cles. The reduced Richardson Hamiltonian should also be an 
appropriate model to describe a mesoscopic paired Fermi sys- 
tem in a trap, with n corresponding to the quantum levels of 
the harmonic trap. Note that in the one-dimensional case, the 
paired states related to each other by time-reversal are identi- 
cal \n) = \tit), while in a three-dimensional potential they are 
\n, I, m) and \n, I, —m) = T\n, I, m). The exact wave-function of 
reduced Hamiltonian ( fTTb is 



|RBCS)exact - C' N J~~[ 



If - F 



|vac), 



where hi = a\cfi 



. . , , and E v are solutions of a system of cou- 

«T "tI y 1 — 1 1 — I J 

pled Richardson's equations Hi 21 1 1 3L 1 1411 . The latter can not 



be solved analytically, but one can show that in the limit of 
large but finite number of particles and small level spacing the 
exact wave-function is indistinguishable from the variational 
BCS wave-function projected to the fixed number of particles. 
This projected BCS wave-function has the form 

|PBCS> = Ca, f ^e-^ N Y\(u n + e^v n ht)\yac), (12) 

" n 

where Cm is a normalization constant and v n and u„ are varia- 
tional parameters. We now use this PBCS wave-function (fT2l 
to calculate the dynamics of the density, which follows after 
the trap and the BCS interaction are turned off. For simplicity, 
we study the one-dimensional case [generalization to higher 
dimensions is formally straightforward; see also Eqs. (|6j [8j 
and the corresponding discussion]. From Eqs. (f2| and ( fT2b . 
we obtain after some algebra 



<n(x, 0> 



2ZAW*,u(r)]r 2 n f m 

Z n In 

{n) N n£{n) N 



(13) 



where {n} N implies all possible sets of N quantum levels, 
{m} n N , denotes all possible sets of N- 1 quantum levels except 
n, f„ = v 2 /u 2 , and if/„[x, u>(t)] is the oscillator wave-function 
with a rescaled frequency cj{t) = w/(l + a) 2 t 2 ). In the limit 
w « A, the sum £ f\ f m , which runs over all possi- 



m " , me m 



depend on the particular level, which is eliminated from it. 
Therefore in Eq. ( fT3l , it can be treated as a constant and deter- 
mined from the self-consistency equation J dx{n(x, f)) = 2N, 
which leads to the following result 



(n(x, f)> = 2N 



00 
n=0 



(14) 



We see that in contrast to the case of a Fermi gas, the sum does 
not have a sharp cut-off at the Fermi level. A smooth cut-off 
is provided by the factor f„ = v 2 /u 2 . E.g., if we assume the 
following expressions for the variational parameters 



1 - 



fn = 



1 + 



(15) 



the sum will be regularized at large indices n » A/a>, as 
f„ ~ I In 2 [we assume here that A is of the order of the 
bulk superconducting gap, A ~ N^ d u>exp(-1/ A), where d is 
the dimensionality]. Note that the only qualitative difference 
between the time-dependent density profile of an expanding 
Fermi gas (0 and the BCS state ( TT4T > is the absence of oscilla- 
tory features in the latter. However, we should note that these 
Friedel oscillations are also absent if the expanding phase is 
initially a strongly interacting but unpaired Fermi liquid. In- 
deed in this case, the initial state is a Fermi sea of Landau 
quasiparticles, which are not the original fermions. Therefore, 
if the typical interaction energy (or temperature) is larger or of 
order u, the oscillatory features are smeared out. Therefore, 
the time-dependent density of a Fermi liquid state and a paired 
state are qualitatively indistinguishable. 

To derive the correlation function, we need to calculate 
the following average <PBCS| a^a^a^a^ |PBCS>. We 
note here that the anomalous pair correlators trivially van- 
ish (PBCS| fl/ ir fl no -' |PBCS) = 0, because the operators in the 
pair correlator do not conserve the number of particles while 
the PBCS wave-function does so by construction. This how- 
ever does not imply that the quartic average has no anomalous 
terms. This is due to the fact that the quartic and higher-order 
averages can not be decomposed into pairs using Wick's the- 
orem if the underlying Hamiltonian is interacting. However, 
the quartic average can be calculated directly using Eq. fT2l) . 

At this point, we introduce the following notation 

a n (x, t) = e'^Vn [x, o)(t)] , where tan (p{t) = a>t. (16) 

This function differs from the usual oscillator wave-function 
only by a time-dependant phase-factor. Using this notation, 
we can write in the limit u> <K A, 



(n(xi)n(x 2 )} = 6(x l -x 2 )(n(x l )}+2- 



-w-i 



1/2 



ble (N - 1) single-particle levels except \ri), does not strongly 



Z N -2Zn . . ... 2Z W _ 2 
+ —j <«(*l)><«(*2)> 

2Z A ,_ 1 Z N 



(17) 
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where Z# = 2 Yl fn- The second term in Eq. (fTTT i repre- 

sents the "anomalous" component of the density, which leads 
to the anomalous spatial correlations in the large-f asymptotic 
behavior discussed previously by Altman et al. |U. Let us 
study the dynamics of this anomalous correlator in detail 



(n(xi,t)n(x2,i)) U 



= 2- 



. Zat-i 



2^ a„(x u t)a„(x 2 , i)fl p - 



(18) 

At f = 0, the phase factor in the function a n (x, 0) JToT l is 
equal to zero 0(0) = 0. If there were no factor J% in the 
sum in Eq. (fT8l . the sum would reduce to Yin 4 r n(x\)^ n (x2) = 
6{x\ - x 2 ), which is a same-point delta-function via the reso- 
lution of unity. The factor which contains information 
about pairing, regularizes the sum at large indices and there- 
fore smears out the delta-function leading to a peak of finite 
width. This has clear physical explanation: At t — 0, the sys- 
tem is in the PBCS paired state with fermions \n, f) and \n, I) 
paired up. The considered limit cj <K A implies that the size of 
a Cooper pair (coherence length, £cp ~ »f/A ~ yjtoN /mA 1 ) 
is much smaller than the size of the trapping potential at the 
Fermi level £ ~ y/N/mcj. In our model, the coherence length 
is small £ » £cp but finite and this leads to a finite width of 
the peak in the correlation function. At larger times t > 0, 
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FIG. 1: (Color online) The figure shows the anomalous density- 
density correlator [see Eq. <18H in rescaled units: C anom (xi, X2', t) = 
a- 2 (0<«[jcicr(0, t]n[x 2 cr(t), f]>, where <r(t) = Vl + co 2 t 2 . The cor- 
relator is plotted as a function of the coordinate x 2 with X\ fixed. 
Note that at t = 0, the correlator is peaked at X2 = Xi, which re- 
flects the fact that the size of a Cooper pair is much smaller than the 
typical lengthscale determined by the trap. At intermediate times, 
the correlator shows a smaller peak, which drifts toward the point 
X2 = -X] . At large times, the rescaled correlator shows a sharp peak 
at xi = -jci, which corresponds to the initial anti-correlation in mo- 
mentum space. 

we should consider the sum of type Ym x l J n{x\HlJn{x 2 )e' m ^ t \ 
which no longer has the form of a resolution of unity. The 
behavior of the anomalous correlator (fT8l was investigated 
numerically (see Fig 1.) using the Ansatz dl5t . We note 
that Fig. 1 shows the evolution of the anomalous peak in 
terms of rescaled variables as defined in the caption. From 



Fig. 1, we see that after the expansion < ut «: A/w, the 
initial peak at x\ ~ x 2 becomes smaller and moves toward 
X\ ~ -x%. At large times t » A/a> 2 , the peak re-appears 
exactly at x\ = -x 2 in accordance with Ref. yl. In our 
formalism, the origin of the peak is due to the fact that the 
phase factor in Eq. (IT6b approaches <p(t — > 00) = n/2 at large 
times. Therefore the sum in the anomalous correlator reduces 
to 2„(-l)"tA„(.x:i)i/',,(x2) ^[fn. Using the symmetry properties 
of the Hermite polynomials, H„(-x) = {-l)"H„(x), one can 
see that in the absence of the factor the sum reduces to 
a delta-function S(x\ + x 2 ) at the opposite points. The "BCS 
factor" smears out the delta-function singularity; again 
due to a finite size of a Cooper pair in the original condensate. 

To summarize, we have derived analytic expressions for 
the dynamics of particle density and density-density correla- 
tion function following a release of a Fermi system from a 
harmonic trap. We considered two types of initial states: A 
Fermi gas and a paired BCS state. Our results for the paired 
state are valid only in the limit of small level spacing co <K A. 
The opposite limit (in which the Cooper pair size is larger 
than the typical length-scale of the trapping potential) should 
be similar to the case of an ultra-small superconducting grain 
(see B12I1 and references therein). In the latter case, there is 
no true superfluidity, but there exist a variety of interesting 
mesoscopic correlation effects (see, e.g., Ref. [17]), which are 
in principle accessible in small atomic systems. 
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